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Abstract 

The fully frustrated XY model with Villain interaction on a square 
lattice is studied by means of Monte Carlo simulations. On the basis 
of the universal jump condition it is argued that there are two distinct 
transitions in the model, corresponding to the loss of XY order and 
Z2 order, respectively. The Kosterlitz-Thouless (KT) transition is 
analyzed by finite size scaling of the helicity modulus at lattices of 
size L = 32 through 128, giving Tkt ~ 0.8108(1). 

The vorticity-vorticity correlation function is used to determine 
two different characteristic lengths, the Z2 correlation length ^, and 
the screening length A, associated with the KT transition and free 
vortices. The temperature dependence of ^ is examined in order to 
determine Tc and the correlation length exponent, i/. The exponent 
is found to be consistent with the 2D Ising value, z/ = 1, and the 
obtained critical temperature is Tc = 0.8225(5). The determinations 
of both ^ and u are done carefully, first applying the techniques to the 
2D Ising model, which serves as a convenient testing ground. 
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1 Introduction 



The critical behavior of the fully frustrated XY (FFXy) model has received 
much attention during the last decade. This is due to the two kinds of sym- 
metries present in the systems and the associated possibility of new critical 
behavior. But in spite of the large number of papers, there is still no con- 
sensus about the phase transition(s) in this model. 

The ordinary frustrated XY model (with cosines interaction) is governed 
by the Hamiltonian 

H = -jJ2cos{e,-ej-Aij), 

{ij) 

where i and j enumerate the lattice sites, 9i is the angle at lattice point i, Aij 
is the quenched vector potential, and the sum is over nearest neighbors. The 
frustration is determined by the sum of Aij around a plaquette (see below). 
In the fully frustrated case - examined in the present paper - this sum is 
equal to vr. 

The peculiarities of the fully frustrated models stem from the two dif- 
ferent symmetries. Beside the rotational symmetry of the XY model, the 
model also has an Z2 symmetry associated with the chirality. The square- 
lattice version in the ground state has a checkerboard pattern of plaquettes 
with positive or negative chirality, corresponding to clockwise or counter- 
clockwise rotation of the spins|jl|]. This is the same symmetry as in the 
anti-ferromagnetic Ising model. 

Some other realizations of YYXY models are the anti-ferromagnetic XY 
model on a triangular lattice]^ , the Coulomb gas with half-integer charges^], 
and the 19- vertex version of the YYXY model[Q. All these models are gen- 
erally assumed to have similar critical behavior. This is also expected to be 
true for the XY Ising model for certain choices of some parameters. 

In the first MC simulations!^ Teitel and Jayaprakash found a divergence 
in the specific heat, consistent with an Ising transition, accompanied by a 
steep drop in the helicity modulus. But since the data did not allow for any 
precise determination of the critical temperature(s), the authors suggested 
two possible scenarios: 

1. At temperatures closely below Tc, the Ising excitations give rise to a 
steep drop in T. When T approaches 2T/7r, the vortex excitations take 
over and produce an universal jump at a temperature Tkt < T^. This 
means two distinct transitions in well-known universality classes. 
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2. The Ising excitations give a transition with an associated non- universal 
jump at the same temperature as the peak of the specific heat. This 
alternative is a single transition in a new universality class. 

While these often have been considered the main options, several other pos- 
sibilities have also been suggested in the literature, as e.g. an Ising transition 
at a lower temperature than the Kosterlitz-Thouless (KT) transition. 

Over the years there have appeared several reports of MC studies where 
the losses of both XY and Z2 orders have been studied. For the anti- 
ferromagnetic XY model on a triangular lattice 0, a study of the heat ca- 
pacity and the XY susceptibility suggested two distinct transitions, with 
Tkt < Tc. The temperature difference was, however, quite small and the 
possibility of a single transition could not be ruled out. In a second study 
of the same model including somewhat larger lattices and with a careful 
analysis of the Z2 transition, the transitions were found to be even closer 
together. The temperature difference was in this study well below the sta- 
tistical uncertainty. The results were therefore suggestive of a single critical 
point. 

Likewise, there are conflicting results for the Coulomb gas with half- 
integer charges in the literature. The flrst study, by Thijssen and Knops[^, 
suggested coinciding transitions whereas Grest found two distinct transitions H 
The conflicting values of the temperature for the loss of XY order was ap- 
parently due to different methods to locate the transition. In the flrst case 
Tkt was determined from the maximum flnite size dependence in 1/e, and 
in the second case from the crossing of 1/e for different system sizes. The 
latter method gives a lower value of the transition temperature. The different 
results for the Z2 transition seem to be due to differences in the MC data. 
Whereas the earlier study reported a drift in position of the peak in the heat 
capacity to lower temperatures with increasing lattice size, such a size depen- 
dence was not verifled in the latter simulation. This discrepancy gave higher 
values of T^. Similar conclusions were also obtained from recent simulations 
of the Coulomb gas with half-integer charges [llO|. The results from this study 



were two distinct transitions; at a lower temperature a KT transition with 
a non-universal jump, followed by a Z2 transition with non-Ising exponents 
(see below). 



As discussed in Ref. [jlT| the dielectric constant at smallest possible wave 
vector k = 2tt/L from the CG simulations, is not an ideal quantity for lo- 
cating the KT transition. With this kind of boundary conditions, there are 
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two finite-size effects working in the opposite directions. That this quantity 
is more or less size-independent only means that these two effects happen 
to nearly cancel each other. This casts doubt on both the KT temperatures 
and the non- universal jumps found in the above studies. 

In order to circumvent the difficulties associated with a precise determina- 
tion of Tkt it has been argued that a determination of the critical exponents 
for the Z2 transition by means of finite size scaling, would be the best way 
to arrive at some firm conclusions. This, at first, seems as a good idea since 
the study of finite size effects right at usually is the by far most efficient 
way to extract the critical behavior by MC simulations. 

In this spirit the correlation length exponent has been determined in 
a fairly large number of studies by means of finite size scaling at T^. In 
the MC simulations |]12|, |13|, this exponent is extracted from the tem- 
perature dependence of various kinds of measures of the distribution of the 
staggered magnetization. The same exponent has also been obtained from 
transfer matrix calculations [jl5|. The results are generally in favor of non- 
Ising exponents, v = 0.85(3), 0.813(5), 0.875(35), 0.80(4), and the critical 
temperatures TJJ = 0.455(2), 0.454(2), 0.4206, 0.454(4). Determinations 
of the same exponent in the 19- vertex version of the FFXy model [Q, the 



Coulomb gas with half-integer charges |[10|], and the XY Ising model[g, \L6 
gave u = 0.77(3), 0.84(3), 0.85(3), and 0.79, respectively. It seems, how- 
ever, to be the case that such finite size scalings in many cases are not quite 
satisfactory and therefore not conclusively, |T3], PH] . 

With the steadily increasing computational resources it has been possible 
to obtain data with high precision for increasingly larger lattices. A recent 
paper reported results for the helicity modulus at a L = 128 system for the 
first time[0. These data has far-reaching implications since it was shown 
that the helicity modulus crosses the universal line, 2T/7r, at a surprisingly 
low temperature, well below the temperatures quoted above for the Z2 tran- 
sition. This must be considered very strong evidence that the XY order 
is lost at a temperature below the Tc obtained from finite size scaling; and 
thus exclude the single transition scenario. However, with this position the 
non-Ising exponents become problematic. The non-Ising exponents are usu- 
ally explained as an effect of the interaction between XY and Ising critical 
excitations - a reasonable explanation only if the two kinds of order are lost 
at the same temperature. 

We therefore have two pieces of evidence pointing in opposite directions. 
The presence of non-Ising exponents strongly suggests a single transition. 
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whereas the early drop in the hehcity modulus seems to exclude this possi- 
bility. 



A consistent view of these matters was recently suggested in Ref. |T7 
The key observation is that a consequence of a KT transition below Tc would 
be the presence of a finite but large screening length A, at the Z2 transition 
temperature, Tc. For finite size scaling to be valid, it is necessary that L be 
much larger than all other finite length scales in the system, and in particular 
A. The large value of A could therefore invalidate earlier finite size scaling 
analyses. The condition for a successful application of finite size scaling at 
Tc, L ^ \, may imply very large systems. 

In the most ambitious study so far of the XY Ising model by means of 
Monte Carlo transfer matrix calculations on infinite strips with widths up 
to 30 lattice spacings. Nightingale et al. again found evidence for non-Ising 
exponents||16|. They did, however, also find an 'internal inconsistency' in 



two different determinations of the thermal exponent ut, which led them 
to call in question the applicability of scaling theory. This inconsistency is 
certainly in line with the suggested failure of finite size scaling due to the 
finite screening length A. 

The main results in Ref. were a precise determination of Tkt, together 
with a demonstration that the staggered magnetization is, indeed, influenced 
by the screening length A, unless L ^ X, i.e. the helicity modulus T ^ 0. In 
order to show that the behavior is consistent with the Ising exponent 1^ = 1, 
the behavior of the correlation length was also examined. This part of the 
study was, however, hampered by two different complications. First, the cor- 
relation function did only fit nicely to an exponential decay for temperatures 
pretty far away from Tc. Second, it was difficult to include the effect of the 
spin waves in an entirely convincing manner. While it is certainly possible 
to argue in favor of the employed technique [|T7[, this is at best only an ap- 
proximative way to compensate for the temperature-dependent effects of the 
spin waves. 

One of the aims of the present study is to improve on the problematic 
points in the temperature dependence of the correlation length. The compli- 
cations with temperature-dependent effects of the spin waves is taken care 
of by performing simulations in an FFXY model with Villain interaction - 
the model dual to the CG with half-integer charges. The point is that both 
the vortex interaction and the vorticity (±1/2) are manifestly temperature- 
independent in that model. To find a reliable technique for determinations of 
the correlation length, we compare with the behavior in the 2D Ising model. 
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In that case we benefit from tlie dual advantages of a fast cluster algorithm 
and exact knowledge of the critical behavior. From the simulations of the 
2D Ising model we show that the region where the true Ising exponent may 
be found is very narrow. A similar analysis of the FFXY model gives the 
same kind of conclusion, and it therefore seems that the data points in Fig. 
5 of Ref. [llTI actually are outside the critical region. 



The main result of the present analysis of the fully frustrated 2D XY 
model with Villain interaction is the existence of two distinct transitions. An 
ordinary Kosterlitz-Thouless transition at Tkt/J = 0.8108(1) followed by an 
Ising transition at Tc/J = 0.8225(5), about 1.4 % above. This is a fairly 
small temperature difference but, as we will see below, the conclusion of two 
distinct transitions is not built on an estimate of the temperature difference 
between two separate transitions. Section p.3.2| gives a strong argument for 
the existence of two transitions which is not based on the determinations of 
the two transition temperatures. 

The organization of the present paper is as follows: In Sec. II we define 
the model, describe the quantities measured in the simulations and some of 
the analyses to be performed on the data. Section III begins our analyses 
of MC data. We shortly describe the MC procedure employed to obtain the 
data and some checks used to validate the results. The major part of Section 
III gives the results from various analyses that take advantage of the finite 
size dependence in the MC data. Among these are the new argument for two 
distinct transitions, the determination of Tkt through finite size scaling of 
the helicity modulus, and an analysis of Binder's cumulant for the staggered 
magnetization. 

Section IV contains the determinations of the characteristic lengths ^ and 
A from the correlation function. In this paper ^ denotes the correlation length 
associated with the Ising-like degrees of freedom, whereas A is the screening 
length associated with the KT transition, which (besides a constant factor) 
is equivalent to the XY correlation length. Since the finite size effects in 
this context are unwanted complications, we take some pains to examine the 
appearance of finite size effects. In order to test some techniques for the 
analysis of correlation functions, and the critical behavior from the correla- 
tion length, we make use of the 2D Ising model as a testing ground. After 
these preliminaries we employ these techniques to the correlation functions 
from the FFXY model to determine the temperature dependence of both ^ 
and A above T^, and C, at low temperatures. This section also contains an 
examination of the effect of domain walls on the vortex interaction. 
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Finally, in Section V we put our results in relation to some results by Berge 



et al. for a model with a variable coupling for one link per plaquette [p^] , and 
summarize our findings. 
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2 Background 



In this section we describe the model, discuss some quantities measured in 
the MC simulations and their relation to the more convenient Coulomb gas 
quantities, and shortly describe some analyses to be applied to the MC data. 



2.1 Model 

The model is defined through the partition function 

r2- dOi 



27r 

where the Hamiltonian for a frustrated system is given by 

(ij) 

In the present case - the Villain version of the FFXY model - the spin 
interaction U{(p), is given by 

oo 



e 



where (f) is an angular difference between nearest neighbors. In the Hamil- 
tonian above Aij is the vector potential, and the frustration is given by the 
rotation of Aij: 

/r = X Ar = — (^r+x/2 " ^r-x/2 " ^r+y/2 + ^r-y/2 

Full frustration, / = 1/2, may e.g. be obtained by setting — everywhere 

and A^ = n at every second row and zero otherwise. Here we introduce the 
discrete difference operator, D = {Dx,Dy), D^fy. = fr+fi/2 — fr-fi/2, and A^ 
and A^ for the vector potential at links in the x and y directions, respectively. 
Associated with the discrete difference is k which is obtained from 

ik^e^^ '' = D^e^^ "" = 2 sin y , 

and also gives = 4 — 2 cos k^ — 2 cos ky. 
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2.2 Measured quantities 



We now describe some of the quantities which are measured in the simula- 
tions, and of central importance for the analyses in Sees. |^ and 



2.2.1 Helicity modulus 

The helicity modulus T, is a measure of the quasi-long range order in XY 
models. It is defined from the increase in free energy due to a small twist A 
across the system in one direction, 

^ d^F 
T = . 

Written in this way, and with current = dF/dA, the helicity modulus may 
be interpreted as the proportionality constant between the applied twist and 
the obtained macroscopic current. 

In MC simulations the helicity modulus is obtained from the correlation 



function 19 




T = 

where Jo = (f/"(0)), and the sum in the second term is over all links in one 
direction, here the x direction. 



2.2.2 Vorticity 

Beside the helicity modulus the main quantity measured in our simulations 
is the Fourier transform of the vorticity. The vorticity is defined in terms of 
the rotation of the current f/'(0jj) = U'{9i — 9j — Aij) around a plaquette 0, 

V = ^iU'{<Pl2) + U'{<j)23) + f/'(034) + f/'(04l)). (l) 

The factor 27r in the denominator is chosen to give the zero-temperature limit 
V = ±1/2. This follows from the angular difference = ±7r/4 in the ground 
state. The steps in the simulations consist of measuring Vr at each plaquette, 
Fourier transforming, 

r 

and accumulating the Fourier components squared, |fkP- 
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It is also common to define the vorticity in terms of the angular differ- 
ences. That corresponds to the chirality, in the context of FFXy models. 
However, an appealing feature of the vorticity defined in Eq. (|Ij) is that it is 
related to some derivatives of the free energy. This is also the reason for the 
existence of some exact relations between the measured vortex correlations 
and the correlations in the 2D CG with half-integer charges. The chirality, 
on the other hand, is somewhat peculiar in that it jumps discontinuously as 
a function of the angular differences. 



2.2.3 The staggered magnetization 



For the study of the Z2 degrees of freedom associated with the symmetry of 
the anti-ferromagnetic Ising model, a convenient quantity is the staggered 
magnetization 

2 / 



M = — 



L2 



(2) 



where the sum is over all the plaquettes of the system, and the alternating 
sign is include to take care of the checkerboard pattern. The factor of 2/1^^ 
is chosen to give M = 1 in a well-ordered system. In an infinite system this 
quantity has a finite value in the low-temperature phase and goes to zero as Tc 
is approached from below, but this sharp behavior is considerably smoothed 



in the finite systems of the MC simulations. Since (— 1 



is directly related to the k = (tt, tt) component of the vorticity. Also useful 
are some powers of the staggered magnetization. 



MP 




Binder's cumulant in Sec. 2.5 is defined from (M ) and (M ). 



2.3 Duality relation and the correlation function 

It has been argued that both the vortex interaction and the average vor- 
ticity at a plaquette are temperature-dependent in the YYXY model with 
cosines interaction |]T^. To avoid this kind of complicating factors one would 
rather have results from the CG with half-integer charges, since both the 
average vorticity and the vortex interaction in that model are manifestly 
temperature-independent. However, since simulations of that model are con- 
siderably more time-consuming, we instead perform simulations of the spin 
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model with Villain interaction, and make use of an exact relation between 
the measured vorticity correlations and the corresponding correlations for 
the CG half-integer charges. In this section we shortly discuss the duality 
transformation, define the Z2 correlation function, g{r) and g(k), and derive 
the link between our measured vorticity correlations and g{k.). 

In the Appendix we discuss the duality transformation pO|, 
the FFXy model. This gives the Hamiltonian 

i/CG _ _4^2 jl j2 mrG{r - rVr', (3) 

r,r' 

where rrir are half-integer charges, rriy. = ±1/2, ±3/2, . . . , and G{r) - the 
lattices Green's function - is the solution to D ■ G{r) = 6r, with G(0) = 
and proper boundary conditions. An excellent approximation for r > 1 is 
27TG{r) = Inr + const. 

It is now convenient to define the Z2 correlation function in terms of the 
CG charges rrir- For the correlation function in ordinary space we write, 

^(r) = 4(-l)"-+''^ (momr) , (4) 

where the prefactors, again, are for normalization and the checkerboard pat- 
tern in the well-ordered ground state. In the low-temperature phase this 
quantity has a finite value in the r ^ 00 limit, whereas it approaches zero 
above Tc. The approaches to these limits are exponential, governed by the 
correlation length ^. As discussed below, the correlation length is, however, 
better determined from the Fourier components g(k). 
The Fourier expansion of the correlation function is 

where we introduce q = (vr, vr) — k and g{q). Together with Eq. (||) this gives 

4 

c/(k) = j-^ (mkm_k) . (5) 

A link between the vorticity correlation function and the correspond- 
ing correlation function of CG charges may be obtained by considering two 
different expressions for the wave- vector dependent helicity modulus T(k), 



21[ applied to 
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which is equivalent to the dielectric function J/e(k)|2^. In the Coulomb gas 
picture, with the interaction 471^7(7 (r) ~ 27rJlnr, T(k) becomes 

T(k) = J - (mkm^k) , (6) 

whereas the same quantity in the XY variables is 

T(k) = Jo - {v^v.^) . (7) 

In the Appendix we show that these two expressions are the same derivative 
of the free energy, which means that they have to be equal. This gives the 
desired link between the correlation functions for our measured vorticity and 
the half-integer variables of the CG model, 

, , , T{J-Jo)L^k^ 
(mkm_k) = [vi^v-i,) + — (8) 

Together with Eq. this gives the desired expression for the correlation 
function g(k) in terms of the measured correlations (fk'W-k)- This is the 
procedure used to determine ^'(k) which is analyzed in Sec. ^. 

The last term in Eq. (||) is under certain conditions very small beside the 
vorticity correlation term (fkf-k)- This is especially the case at q ~ for 
temperatures around T^, which means that the determinations of ^ and i' in 
Sec. ^ would be influenced only very slightly by neglecting this correction. 

The relation between the CG correlations and the vorticity correlations 
in ordinary space has been discussed in Ref. [^. For the case with Villain 
spin interaction the result was (momr) = (voVr)- However, as seen in their 
derivation this holds for |r| > 1, only. For general r the Fourier transform of 
Eq. (IP gives 

4T(J- Jo)/(27rJ)2, ifr = 0, 
(momr) - (voVr) = { -T{J - Jo)/(27rJ)2, if |r| = 1, (9) 

0, otherwise. 



In Sec. p.2| the above relation for r = is used in a consistency check. 



2.4 Analysis of T 

In this section we discuss the size-dependence of the helicity modulus T and 
its relation to the universal jump and Kosterlitz' renormalization group (RG) 
equations. We will focus on the dimensionless quantity T/T. 
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For the finite-size scaling analysis of T/T one assumes that the size- 
dependence in this quantity is related to the behavior of a set of RG trajectories J 
The starting point in parameter space, and thereby the relevant trajectory, is 
determined by the temperature. These trajectories behave differently in the 
low- and high-temperature phases. In the low-temperature phase they termi- 
nate at finite values of T/T whereas they continue to infinity, corresponding 
to T ^ 0, in the high-temperature phase. 

The last trajectory in the low-temperature phase ends at the universal 
value T^G^Y ^ T/{2nT) = 1/4. This means that the helicity modulus for 
that very temperature in an infinite system is T = 2T/tt. The jump of this 
quantity to zero is the well-known universal jump||2^, |25|| . 

The abrupt universal jump of an infinite system is, of course, not seen in 
finite systems. Since T decreases with increasing system size, the universal 
jump conditions T = 2T /ti may, however, be used to establish an upper limit 
for Tkt- The temperature obtained in that way is a rigorous upper limit, 
since the universal jump condition constitutes an absolute stability criterion. 

The approach to the universal value, T7r/(2T) = 1, with increasing sys- 
tem size, may also be used to examine the critical properties. From Koster- 
litz' RG equations the finite-size scaling relation for becomes [p^ 



Kosterlitz' RG equations are expected to be valid only in the limit of low 
vortex density. This means that the above finite size scaling relation is ex- 
pected to be valid only at low renormalized vortex density. Accordingly, 
Ti7r/(2T) should be not too far from unity - T renormalized out to length 
scale L should not be too far from the fully renormalized T out to infinity. 
This implies the dilute limit for sizes bigger than L. 

The same idea may also be used both above and below Tkt- A more 
complete discussion is given in Ref. Close to Tkt we expect 



^ = l + ccoth[2c(£o + lnT)], T < T^, (11) 
^ = l-ccot[2c(4-lnT)], T>T,., (12) 



where (.q and c are free parameters to be determined from the fits. Eq. (10) 



is the c — >■ limit of Eq. (|TTI). c vanishes as Tkt is approached from below 
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or above as [OB 



c = 5^|T/Tkt-1|. 

In the high-temperature phase io is identified with the logarithm of the 
screening length, A. In the immediate vicinity of Tkt the temperature de- 
pendence of io should therefore be given by Kosterlitz' expression p3|] 

C , , 

= + const, (13) 



/T/Tkt 

where [E2 



2.5 Binder's cumulant 

Binders' cumulant is a convenient quantity that, in most cases, facilitates 
determinations of both the critical temperature and the correlation length 
exponent v. Even though the quantity was originally presented in terms of 
averages over blocks of different size in a single simulation with a fixed total 
system size[^, it may also be used with data from systems of different size. 



Binder's cumulant is obtained from some moments of the order parameter, 

(15) 

The crucial property of U is its size-independence precisely at T^. Therefore, 
plotting U versus temperature for several different sizes is expected to give 
an unique crossing point at the critical temperature. Furthermore, the cor- 
relation length exponent v may be determined by plotting the data against 
(T — Tf)L^I^ . The correct value of v is expected to give a collapse of that 
data onto a single curve. In practice there are, however, often corrections to 
scaling which make the conclusions from this kind of analysis less direct and 
precise. 



2.6 Boundary conditions 

It has recently been pointed out that periodic boundary conditions (PBC's) 
in the XY model may be generalized by including twist fluctuations along the 
X and y directions in the system - fluctuating boundary conditions (FBC's) 
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TT]. There are several advantages with considering such a generahzation. 
First, it is with these boundary conditions that the Villain version of the XY 
model is exactly dual to the CG with periodic boundary conditions. Second, 
the finite size effects in several quantities work in the opposite way after the 
inclusion of these twists. Finally, with a self-consistently chosen amplitude 
of these twists, the finite size effect on the correlation function turns out to 
be virtually eliminated. 

The self-consistent boundary conditions do not seem to be applicable in 
the fully frustrated case. This is possibly an effect of the Z2 fluctuations. 
It is, however, possible to obtain the correlation function as in an inflnite 
system by taking the average of data for PBC's and FBC's, cf. Fig. 1 in 
Ref. [^. This technique works up to, and possibly slightly above Tkt- 
At higher temperatures both sets of data go down with increasing lattice 
size, which makes it considerably more difficult to extract any result for the 
thermodynamic limit. 
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3 Finite size analyses 



The different methods to analyze MC data may, generally speaking, be di- 
vided into two classes. The most obvious one is to calculate the correlation 
functions and determine the correlation length and the associated exponents 
from this kind of data. In this kind of analyses one is interested in the behav- 
ior of an infinite system and, accordingly, the finite-size effects are undesired 
complications. This kind of analyses is the subject of Sec. 

The second class of methods instead take advantage of the finite-size 
dependence in the MC data. This is generally a more efficient approach to 
analyzing the critical behavior. In this Section we employ some techniques 
that make use of the finite-size dependence in various ways. 

After a short description of the simulations and some checks employed to 
validate the results, we focus on results from the universal jump condition 



in Sec. p.3| . In Sec. |3.4| we perform finite-size scaling analyses of the helicity 
modulus T both right at Tkt and in the immediate neighborhood around 
Tkt- With this determined value for Tkt we then take a closer look at the 
data from the universal jump condition in Sec. |3.5| . 

To obtain a reference temperature we then apply finite-size scaling anal- 
ysis of Binder's cumulant at Tc. Just as in the related models this kind of 
analysis gives u < 1. As suggested in Ref. [|1^ this seems to be an artifact 
of the presence of a finite screening length A associated with the nearby KT 
transition. 



3.1 Monte Carlo simulations 

The Monte Carlo simulations were performed with the ordinary Metropolis 
algorithm with sequantial sweeps over the lattices. One such sweep with one 
trial update per spin is called a MC step. For most of the data there were four 
MC steps between consecutive measurements. But since it was noted that a 
major part of the computer time, especially on the large lattices, was used 
in the Fast Fourier Transform of the measured vorticity v^, the simulations 
for L = 128 and 256 for determinations of ^ were performed with as much 
as 64 MC steps between consecutive measurements. For the latter data the 
number of MC steps are given in Table ^. 

In Sec. ^ we make use of MC data from the 2D Ising model as a con- 
venient testing ground for the methods used to analyze the FFXY model. 
These simulations are performed with Wolff's cluster algorithmpU|. All the 
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simulations were done on a set of DEC-alpha workstations. 



3.2 Monte Carlo data 

A MC study is, of course, never more reliable than the underlying data. It 
is therefore essential to check that the program, indeed, does provide correct 
data. This may be done either by comparing with previously published 
results or by making use of some consistency tests. 

To the best of our knowledge there is no published MC data to compare 
with for the Villain version of the fully frustrated XY model with ordinary 
PBC's. For the case with FBC's it is, however, possible to compare with data 
for the half-integer CG[]TD|. For L = 8 and T/J= 0.82 our simulations give 



Binder's cumulant, U = 0.5786. This temperature corresponds to T*"^ 
0.1305, and as expected our value for U lies right in-between the values for 
U at T^^ = 0.130 and 0.131 in Fig. 3 of Ref. As a second test we 

compare the values of T{k = 2tt/L) for L = 16 and T/J = 0.82. Again, the 
value obtained from our simulations, 0.556788, is in good agreement with the 
corresponding values in Fig. 4 of Ref. fl^ . 

For the bulk of our data, obtained with ordinary PBC's we have to resort 
to internal consistency tests. One such test is suggested by the analogy with 
the CG with half-integer charges. In that case, the charges m = ±1/2 give 
(m^) = 1/4. Actually, the value 1/4 turns out to be a lower bound since 
the CG also includes charges of non-lowest order, i.e. m = ±3/2. For our 
measured quantity v^, there is no such simple result, but as discussed above 
there is an exact relation between these two quantities, Eq. (H). 

The behavior of both (f ^) (squares) and (m^) from Eq. (circles) is 
shown in Fig. |l|. Whereas (f^) is seen to decrease with temperature we 
find that (m^) indeed is very close to 1/4. More precisely, the results are 
(m2) = 0.250000 at T/J = 0.4, 0.250003 at T/J = 0.77, and 0.250078 at 
T/J = 0.9. This constitutes a confirmation of the correctness of the MC 
data. From the minute deviations of (m^) from 1/4 it is also possible to 
obtain estimates of the fraction of plaquettes with non-lowest order charges 
in an equivalent CG|^. For the temperature interval in Fig. |I| the data 



indicates that this fraction would be from 1.5 x 10 to 39 x 10 
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Figure 1: The average vorticity squared for a single plaquette (f^), together with 
the corresponding Coulomb-gas quantity (m^) as functions of temperature. That 
the latter quantity, obtained through Eq. @ , is close to 1/4 is a consistency test of 
the simulations. The small deviations from 1/4 give information about the density 
of non-lowest charges, m = ±3/2, in the corresponding Coulomb gas. 

3.3 Universal jump condition 

In this section we make use of the universal jump condition to obtain both an 
upper limit of Tkt and a strong evidence in favor of two distinct transitions. 
Since the universal jump condition is an absolute stability requirement, we 
believe the argument of the present section to be especially free of objections. 
Whereas the more precise results in the later sections are obtained on the 
basis of additional assumptions, the direct use of the universal jump condition 
is particularly clean. 

For the following discussion we introduce the size-dependent transition 
temperature T^, as the temperature where the helicity modulus for system 
size L is equal to the universal value, = 2T^/7r. 



3.3.1 Upper limit of Tkt 

The universal jump condition was applied to the fully frustrated XY model 
with cosines interaction in Ref 
temperature 



131 to establish an upper limit for the KT 
From the intersection of the MC data for L = 128 with 
the universal value, these authors found, as discussed in the Introduction, 
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Figure 2: Determinations of the size-dependent KT temperatures T^r^. The 
dashed hue is the universal jump condition. The sohd hues are second order 
polynomials in T/J obtained from fits to the MC data. 



T^r^^^/J = 0.449(1) J, clearly below the values of Tc determined with finite 
size scaling. 

The same approach with the data for the Villain version is shown in Fig. 
^ The upper limit obtained for L = 128 is T^t^'' ^ 0.816J. Figure ^ shows 
the size- dependence of T^. The dashed line is from the analysis in Sec. |H75. 



3.3.2 Two transitions 

We now turn to the argument for two distinct transitions based solely on the 
universal jump condition |]2^, |25[| . 



The starting point is that the staggered magnetization for an infinite sys- 
tem, Moo, vanishes at the Z2 transition temperature Tc. To establish the 
existence of two distinct transitions it is therefore sufficient to examine M^o 
right at Tkt- A non-zero value of Moo(Tkt) would be an unequivocal demon- 
stration of Ising order, which then implies that this order is lost at a higher 
temperature, Tc > Tkt- The determination of Moo(Tkt) at first seems very 
difficult since, beside the usual problem of approaching the thermodynamic 
limit, the value of Tkt has to be known with high precision. The universal 
jump criterion used so far, is only capable of yielding upper limits. 

A way around both these difficulties at the same time is to focus on 
Mi (Tkt), ^^e staggered magnetization at finite lattices at the size-dependent 
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Figure 3: The size-dependence of T^. The cross is Tkt from Sec. 3.4 and the 



dashed hne is from the analysis in Sec. 3.5. 



KT temperatures, and examine the behavior of this quantity as a function 
of system size, L. The point is that the desired quantity Moo (Tkt) is the 
large-L hmit of MiiTj^), and that the staggered magnetization is readily 
determined for each L. 

The results from this analysis is shown in Fig. |[ The figure shows that 
the staggered magnetization at is an increasing function of lattice size. 
The figure gives 0.744 as a lower limit of Moo (Tkt), and a naive extrapolation, 
that neglects the curvature, suggests Moo(Tkt) > 0.768. 

We consider this to be a very strong argument that the Z2 order persists 
at the KT transition temperature. For the opposite to be true, this increas- 
ing trend toward a finite value of Moo (Tkt) should change to a decreasing 
trend down to zero. Even though this possibility could never be ruled out 
from the data for finite systems alone, such a change in trend seems very 
unlikely. Furthermore, the more detailed analysis in the following sections 
yields Moo (Tkt) = 0.783(2), entirely consistent with the increasing trend in 
Fig. H It should also be noted that this line of evidence does not depend on 
the assumption of an universal jump. The argument holds equally well with 
a jump T7r/(2T) = g, g > 1. {g < 1 is excluded by stability.) 

The actual determination of ML{Tj:^) is illustrated in Fig. ||. Ml{tI^) 
is obtained directly from the value of Ml where = 2T/7r (dashed line in 
Fig. D. 
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Figure 4: Evidence for two distinct transitions. The plotted quantity is known 
to approach Moo(Tkt) in the limit L ^ oo. The dashed line shows a naive 



extrapolation, whereas the dot is this quantity from Sec. 4.5 obtained with the 



value of Tkt from Sec. 3.4. 




21 



L 






IL 


-Grange/ J 


16 


1.3932(17) 


-8.545 


-19.9 


0.021 


24 


1.3179(20) 


-11.06 


-143.6 


0.015 


32 


1.2731(12) 


-12.86 


-207.4 


0.01 


48 


1.2235(13) 


-15.75 


-190.8 


0.008 


64 


1.1989(13) 


-17.55 


-352.1 


0.0065 


96 


1.1714(15) 


-20.20 


-439.0 


0.0051 


128 


1.1578(20) 


-21.30 


-304.1 


0.0045 



Table 1: Parameters from fitting MC data for the helicity modulus to Eq. 
([T6|). The data included in the fits are restricted to \T / J— 0.8107| < Tj-ange/ J- 

3.4 Kosterlitz-Thouless transition 

The purpose with this section is first to determine the Kosterlitz-Thouless 
temperature, Tkt, and, second, to examine the behavior closely below and 
above this temperature. The method employed is finite size scaling analysis 
of the helicity modulus as discussed in Sec. The basic idea is that the 
size-dependence of the helicity modulus T, at and in the vicinity of Tkt, niay 
be obtained from Kosterlitz RG equations||2^ as given by Eqs. (0)' (HID ^-^^ 



The analysis of the helicity modulus T, is based on a large amount of 
MC data. In order to make efficient use of the data and get T/, as continu- 
ous functions of T, we first determine the helicity modulus as second order 
polynomials, one for each L, in t = T/J — 0.8107: 

^ = «L + /?Lr + 7Lrl (16) 

These second order expansions are only expected to be valid within rather 
narrow temperature intervals. We therefore only include data in the fits 
for temperatures |r| < Trange/<^, where Trange decreases with increasing L. 
The parameters from this analysis together with the size of the temperature 
intervals are shown in Table B.4|. 



3.4.1 Determination of Tkt 

We now apply the finite size scaling relations for the helicity modulus as 



discussed in Sec. |2^ . Since these relations in this system only are expected 
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Figure 6: The quality of the fit from fitting to Eq. ([lO|). Lmin is the smahest 
size included in the fit. For a good fit one expects x^/DOF ~ 1. The figure indi- 
cates that the fit is not very sucessful when including small lattices, but becomes 
acceptable for L^^-a = 32. 



to be valid for fairly large lattices we first follow the procedure in Ref. ||27 
and perform the analysis with systems of size L = Lmin through 128 and 
various values for Lmin. The errors in the fits are shown in Fig. ^ On the 
basis of this analysis we conclude that L^m = 32 does give a good fit. This 
is the same choice as for the FFXF model with cosines interaction in Ref. 



Figure |^ shows the good fit of the MC data to Eq. ( [I0| ) obtained by 
adjusting Tkt and Iq. The obtained value for the KT temperature is Tkt/</ = 
0.8108(1). We consider the good fit to Eq. (PI1|) to be very strong evidence 
for an ordinary KT transition. 

Note that Tkt/ J = 0.8108 is well below the upper limit T/J = 0.816 
from the universal jump criterion in Fig. 3. It is also slightly lower than 
what a simple linear extrapolation of the four lowest points to 1/L = in 
Fig. ^ would suggest. 

3.4.2 Finite size scaling around Tkt 

We now shortly discuss the critical behavior in the immediate vicinity of Tkt- 
The approach closely follows Ref. ||27|| . 

In fitting our MC data to Eqs. (|TT|) and (|l^) we fix the temperature, 
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calculate TL7r/(2T) from the parameters in Table and adjust c and Iq 
to get the best possible fit. /q is the logarithm of the screening length, A, 
in the high-temperature phase. Below Tkt this quantity has no such direct 
interpretation. The results from this kind of fitting for several temperatures 
around Tkt are shown in Table |^. 

Figures ^ show the temperature dependence of c. Just as in the analysis 
of the ordinary XY model the values of the slopes B obtained from the 
low- and high-temperature data are, within statistical errors, the same. Also 
shown in the figures are the corresponding values of C from Eq. (|T^. The 
slope of £o versus 1 / ^T/T^j, — i in Fig. is in good agreement with the 
values from Figs, pi and our estimate for the slope becomes C = 0.54 ±0.02. 



In Sees. |3.5| and |4.4| we will obtain different values for C, but we consider the 
present determination to be the more reliable one for two reasons. First, it 
does build on excellent agreements with predictions from the Kosterlitz' RG 
equations and, second, in contrast to other determinations, this method does 
probe the behavior in the immediate vicinity of Tkt- 

3.5 Size-dependence of 



In Sec. ^]3| we made use of the universal jump condition to determine a kind 
of size- dependent KT temperatures Tkt^ as upper bounds for Tkt- From Fig. 
Rl it seemed difficult to extrapolate such data to the thermodynamic limit. 
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VT/Tkt - 1 



Figure 8: The temperature dependence of the parameter c in Table ^. The data are 
obtained by fitting Tl to Eqs. (0) and (|T2|) , respectively, (a) Low temperatures, 
T < Tkt- In this case c = Toovr/(2T) — 1, and the figure therefore illustrates the 
approach of Too to the universal value as \/l — T/Tkt- At high temperatures, 
panel (b), c has a similar square-root cusp, though it is no longer related to Too. 
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13.66 


1.13 
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1.66 


0.814 


0.1775 


10.45 


1.98 


0.815 


0.2029 


9.34 
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Table 2: Parameters from the fitting of MC data to Eqs. ([TT|) and ([T^ ) below 
and above Tkt, respectively. 



We now demonstrate that the size-dependence of has the same form 
as the Kosterlitz' expression for the temperature dependence of the correla- 
tion length, i.e. 

C 

InL = const H — , (17) 
/tS/Tkt - 1 



Figure |TO|(a) shows the approach of t'^^ to Tkt = 0.8108J with increasing 
L. The points do, indeed, fall on a straight line. This is the same dashed 
line as in Fig. ^ It should be noted that it not seems possible to link 
this behavior directly to the divergence of the screening length, A. In the 
present approach the slope is C ~ 0.265, whereas the same constant in the 
temperature dependence of /q yielded C ~ 0.54 which is about twice as big. 

This size-dependence of T^^^ is also found in the ordinary 2D XY model 
with no frustration. This is illustrated in Fig. |lO](b). Here the slope is ~ 0.85, 
and, again, the temperature dependence of the characteristic length gives a 
slope that is about twice as big[^ . 
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Figure 9: The temperature dependence of the screening length above TkT) pre- 
dicted by Kosterhtz, Eq. (13). Beside an additive constant, £o is equal to the 



logarithm of the screening length. This is a very direct determination of the slope 
C. The data, which is also listed in Table ^, was obtained by fitting to Eq. 



3.6 Binder's cumulant 

As discussed in the Introduction the existence of two transitions close to each 
other may give rise to problems with ordinary finite size scaling. The evi- 
dence for two transitions given in the previous Section, strongly implies that 
this actually is the case for the FFXY model. This means that the critical 
properties not are accessible with finite size scaling at Tc unless the systems 
employed are considerably larger than the correlation length associated with 
the other (here the KT) transition. 

The purpose with the present finite size scaling analysis is therefore not 
to extract the correct critical behavior, but rather to provide a reference 
temperature and to verify that the FFXy with Villain interaction indeed 
does behave in a way that is similar to the more studied FFXY model with 
cosines interaction. 

Figure |ll|(a) shows Binder's cumulant versus temperature for L = 8, 16, 
32, and 64. As discussed in Sec. pTSj U is expected to be size-independent 
right at the critical temperature. This is not quite borne out by the data. A 
close look reveals that the crossing points move slowly to lower temperatures 
for larger system sizes. For the pairs of lattice sizes L = 8, 16, L = 16, 32, 
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Figure 10: The size dependence of T^^^. (a) For the YYXY model, with lattice 
sizes L = 16, 32, 48, 64, 96, and 128. The value for the transition temperature is 
taken from the determination in Sec. (b) The corresponding quantity in the 
ordinary XY model with no frustration. This is with / J = 0.8921 p7|| . Note 
that the slope C in these cases are entirely different from the corresponding slopes 
from the temperature dependence of the screening length. 
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and L = 32, 64, the crossing temperatures are 0.827, 0.825, and 0.824, re- 
spectively, though the two last temperatures are within the statistical uncer- 
tainties. This is in good agreement with 27r x 0.1315(3) = 0.826(2) obtained 
from simulations of the CG with half-integer charges on lattices with size 
L = 10 - 24 in Ref. |]10|- Since the ordinary periodic boundary conditions in 



that simulation corresponds to FBC's in spin models (cf. Sec. |2.6|) , whereas 
the present simulations are performed with PBC's, the value of the cumulant 



at criticality is, however, not expected to be the same[|T7 



Figure |TT](b) shows the data collapse. Following Ref. we assume 
that Ul{T) = (pitL^/"), where t = T/T^ - 1. We then expand 0(x) = 
00 + 012^ + 023;^ for small x and adjust these three parameters together with u 
and Tc to get the best possible fit. With data close to T^, (0.605 <U < 0.640) 
for system sizes L = 16, 32, and 64, we obtain Tc = 0.8244 and u ~ 0.81. 
This value of u is in good agreement with the published values, listed in the 
Introduction. 
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Figure 11: Binder's cumulant, Ul, for 
less size-independent for L > 16 at T/J 
Just as in tlie FFXY model with cosines 
u<l. 



different lattice sizes, (a) U is more or 
~ 0.824. (b) An attempted data collaps. 
interaction this kind of analysis suggests 
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4 Correlation lengths 



In the previous Section we performed a number of analyses with methods 
that take advantage of the finite size effects in the MC data. This is usually 
the most efficient way to determine the critical behavior from Monte Carlo 
simulations. The alternative approach is to determine the correlation length 
from the length-dependence of some correlation functions and then extract 
the critical behavior from its temperature dependence. In the present section 
we take this alternative route. 

In order to test some techniques for determining the correlation length and 
extracting the correlation length exponent and the critical temperature, we 
first present an analysis of the 2D Ising model. Because of the dual advantage 
of a known critical behavior and a fast cluster update algorithm, this model 
serves as a very convenient testing ground. One minor difference between 
the analysis of the 2D Ising model and the FFXy model is due to the anti 
ferromagnetic ordering. Whereas the critical behavior of the ferromagnetic 
Ising model is obtained in the k ^ limit, the corresponding critical behavior 
in the FFXY model manifests itself at k = (vr, tt). This is taken care of by 
performing all the analyses in Sec. in terms of q = (vr, tt) — k instead of 
k, cf. Sec. O. 



In this section we determine two different characteristic lengths from our 
MC data for (fk^-k)- The reason that it is at all possible to define two dif- 
ferent characteristic lengths is related to the above discussion. Whereas the 
screening length A is determined from the k ^ limit of these correlations, 
the correlation length ^ is determined from the limit g — > 0. 

For the determination of the characteristic lengths one like to have the 
correlation function for an infinite system. It is therefore of great impor- 
tance to know when the undesired effects of the finite lattice size set in. 
Before applying the obtained techniques to the FFXF model we do a careful 
analysis of the finite-size effects in this model. The results corroborate the 
suggestion |]17| that the correlation function is plagued by finite size effects 
unless the system is large enough that T ^ 0. 

After these preliminaries we then turn to determinations of the correla- 
tion length ^, the critical exponent z/, and the critical temperature Tc, in the 
FFXY model. Much as expected from the evidence of two distinct transi- 
tions, the behavior is found to be consistent with an ordinary Ising transition, 
u = 1. However, for this demonstration it turns out to be necessary to ex- 
amine the behavior fairly close to T^, which corresponds to large correlation 
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lengths, ^ > 10. The screening length A associated with the KT transition is 
also determined, and its temperature dependence is found to be entirely dif- 
ferent from the behavior of ^, but in good agreement with Kosterlitz' result, 
Eq. (Pm. 

We finally turn to the behavior of ^ at T < Tc. At these temperatures 
we have no data with T ^ 0, (and if we had, there might also be problem 
fulfilling L ^ C,) and we therefore need some other methods to avoid the 
finite-size effects. The solution is to restrict the analysis to T < Tkt, where 
the different boundary conditions of Sec. |2.6| may be employed. However, the 
temperature dependence of ^ below Tc (and Tkt) does not seem to be useful 
for assessing the critical behavior, possibly an effect of the presence of the 
KT transition between the obtain data and T^. 

4.1 The correlation length in the 2D Ising model 

The Hamiltonian of the Ising model is 

{ij) 

where i and j numerate the lattice points, Si = ±1, and the summation 
is restricted to nearest neighbors. In two dimensions the correlation length 
exponent is z/ = 1, and at a square lattice the critical temperature is known 
to be 

T// J = ^ ^ 2.269. 

ln(v^+l) 

Since = 1, a plot of 1/,^ versus T is expected to yield a rectilinear 
behavior down to T/. However, the verification of this turns out to require 
data fairly close to T/, large correlation lengths, and therefore rather big 
lattices. For most purposes this exercise is pointless, since - beside being 
obtained from the exact solution - the value u = 1 may be verified from MC 
simulations by means of finite size scaling at T/. But since this kind of finite 
size scaling does not seem to work in the FFXF models for the accessible 
lattice sizes, we have to resort to analyses of the correlation functions. With 
that background, analyses of the correlation function for the 2D Ising model 
serves as a help to develop techniques for similar analyses of the FFXY 
model. Beside the benefit of the exactly known critical behavior, the analysis 
of the Ising model is greatly simplified by means of the cluster algorithm pO[ 
that is instrumental in obtaining MC with small statistical errors. 
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4.1.1 Determination of the correlation length 

At first sight the obvious way to determine the correlation length is to ex- 
amine the exponential decrease of the correlation function g{r) down to zero. 
This amounts to adjusting the parameters A and C, to obtain best possible 
fit to the expression, 

g{r) = Ae-'/^. 

At temperatures closely above Tc, the above expression should be modified to 
take correlations across the whole system into account. This is customarily 
done by instead fitting to an expression with the periodicity of the system, 

g{r) = A [e-'/^ + e^^^-'^^/^) . (18) 

It is, however, difficult to obtain reliable values for the correlation length with 
this procedure. The main complication is that the optimum value of ^ does 
depend on the range in r employed for the fit. This is not too surprising, 
since the pure exponential decay only is expected for very small values of 

An alternative determination of ^ by means of g{k), the Fourier compo- 



nents of the correlation function, has been suggested in Ref [33 



«'°) 1. (19) 



2ir\|9(27r/i) 

An advantage with this expression is that no fitting is needed, and that 
the arbitrariness involved in choosing the fitting interval is eliminated. This 
expression may be derived from 

The relation to an exponential decay is obtained since the Fourier transform 
of this function is the Bessel-ii^o function, 

g{r) oc Ko{r/^), 

with the limiting behavior ~ e"'"/'^. 

In Fig. we display some determinations of the correlation length in the 
2D Ising model. The open circles are obtained with Eq. (^8|) whereas the 



33 



0.18 



0.16 - 
0.14 - 
^ 0.12 - 

o.io'- 




0.08-.----" 

II ' 

0.00 0.02 0.04 0.06 0.08 0.10 0.12 

1/^min 

Figure 12: Results for the correlation length in the 2D Ising model. The open 
circles are ^ from fitting g{r) to Eq. ( |l^ ) for rmin < r < 2r^in- The solid squares 
are determinations of ^ from the small- /c limit, Eq. (^). The data obtained in 
that way thus seems to correspond to the r^[^ oo limit. The MC data is for 
L = 128 and, from top to bottom, T/J = 2.45, 2.41, and 2.37. 



solid squares are from Eq. (|T^). The curves are for different temperatures, 
from top to bottom, T/J = 2.45, 2.41, and 2.37. The different values of ^ 
are due to the different ranges of r for the data included in the fit. We use 
g{r) for rmin < r < 2r^i^. The x axis shows l/r^in- We find that as rmin 
increases, the inverse correlation length decreases towards the solid squares 
from Eq. (19). This suggests that the k limit really does probe the 
long-distance limit. On the basis of this comparison we believe that Eq. ([19D 
gives a reliable way to determine the correlation length ^. 

A precise determination of ^ with Eq. (|19|) requires fairly long MC sim- 
ulations. This is the case since g{27i/L) and g{0) measure the amplitude of 
the largest fluctuations in the system, with the correspondingly long decor- 
relation times. A way to reduce the effect of statistical errors is to include 
some more fc-vectors in the analysis. This is motivated by the difficulty to 
obtain good accuracy from Eq. (p!9D on data from the FFXY model at large 
lattices (L = 128, 256). However, using g{k) in a too large interval will affect 
the correlation length ^. Assuming that the exponential r-dependence only 
holds for r ^ ^, or, similarly, that the asymptotic /c-dependence only is valid 
for A; ^ 27r/^ we restrict ourselves to making use of data from wave- vectors 
k < tt/^, only. For small values of ^ this is not very restrictive, and since 
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we are interested in the small-fc limit we impose the additional condition 

k < ^/0J. 

The procedure to determine ^ is then to first fit the data to 



1 



(21) 



(where we also include a second order term in to take care of the curvature 
in the data) and then extract the correlation length through 



In the limiting case with only two wave-vectors, this procedure is readily 
shown to be equivalent to Eq. ([T^). 

4.1.2 Determination of the critical behavior 

Figure |l^ shows the temperature dependence of the correlation length in the 
2D Ising model for three different system sizes, L = 64, 128, and 256. In the 
vicinity of T/ the data reveals some finite-size effects. The correlation length 
becomes smaller in a too small system. The present data seems to suggests 
that the determinations are reliable only if L/^ > 5. 

Also apparent in the figure is a slight curvature in the data. The expected 
linear behavior is found only right above T/. Estimates of T/ may be obtained 
by fitting 



with A and as free parameters. In these analyses we only make use of 
data for L = 128 and 256. Due to the curvature in the data, the critical 
temperature obtained in this way does depend on the temperature interval 
for the fit. For temperatures closely above T/ we only include data points 
with L/C, > 5. This gives L-dependent lower limits for the temperature 
interval. The upper limit of the temperature interval is given by Tmax- The 
fit is then performed for several different values of T^ax- The dependence of 
T/ on Tjnax is shown by open circles in Fig. ^^a). The dashed line is the 
exact value of T/. For large T^ax (large temperature intervals) the analysis 
yields too low estimates of the critical temperature, but with decreasing Tmax 
the estimated T/ increases towards the correct value. 
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Figure 13: Correlation length in the 2D Ising model for L = 64, 128, and 256. 
The solid line is from fitting to Eq. (|2|) for T/ J < 2.33 with data that fulfills 
L/e > 5. 



This linear fit presumes a known value of the correlation length exponent 
V. Since the value of v in the YYXY model is highly disputed it is also of 
interest to perform the fit with i/ as a free parameter. That is done by fitting 

l/i = A{T - T,y, (24) 

with u, A and T^, as free parameters. Again we repeat the analysis for several 
different Tmax- The values of T/ and u as functions of T^ax are shown by 
solid squares in Figs, pi We find that both quantities approach the expected 



values as the temperature interval is reduced. The erratic behavior at low 
values of T^ax is due to statistical errors that become significant in such 
narrow temperature intervals. 



4.2 Finite size effects 

In the previous section we found, much as expected, that it is necessary 
to perform the MC simulations of the 2D Ising model at systems that are 
considerably larger than the correlation length. The criterion was L/^ > 5. 
When finite size effects set in we get larger correlations and, thereby, too small 
correlation length. In this section we demonstrate that in the analysis of the 
FFXY model, this condition has to be supplemented by a second one that 
is related to the additional XY degrees of freedom, and the corresponding 
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Figure 14: Determinations of Tc and u in the 2D Ising model from the correlation 
length data by fitting to Eqs. (p^) and (|2^). The open circles are from fits with the 
exponent keeped fixed, 1^ = 1, whereas the solid squares are obtained with both 
Tc and u as free parameters. The general trend is that both quantities approach 
the exact values, indicated by the dashed lines, as the temperature interval is 
decreased. 
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screening length, A. In terms of the hehcity modulus this condition may be 
written T ^ 0, which is equivalent to L ^ A. 

In order to check where the finite size effects become important for the 
determination of ^, it is convenient to monitor g{q) for q ^ from different 



lattice sizes as functions of temperature. This is done in Fig. |T5|. Panel (a) 
shows g{q = 0) for lattice sizes L = 64 and 128. g{q = 0) for the even larger 
system, L = 256, suffers from large statistical errors and is therefore not 
included in the figure. Comparison between L = 128 and 256, are instead 
performed with g{q = 27r/128), shown in Fig. p!5|(b). 

In both these figures the helicity modulus for the smaller size is included 
as solid lines. Comparing the data for the different sizes we find that they 
start to differ at about the temperature where the helicity modulus becomes 
appreciably different from zero. This implies that determinations of the 
correlation length from g{q) are uncertain if the system is not large enough 
to ensure that T 0. 

The suggested link from the two previous figures between the size- dependence 
in g{q) and the helicity modulus may be understood by examining the wave- 
vector dependent helicity modulus T{k) for different system sizes. Implicit 
in this discussion is the close relation between (yf(k) and T(k) in Eqs. and 
(H). Figure |16| shows T(fc) at T/J = 0.84, which is well above both Tkt and 
T 

At each temperature above Tkt the helicity modulus T = T(fc = 0) van- 
ishes for sufficiently large L. A finite value of T may therefore be considered 
a finite-size effect. The message of Fig. ^ is that if T suffers from finite-size 
effects, then the same is true for all the other components T(fc), as well. 
The converse also appears to be true. For system sizes with T ^ 0, T{k) 
is independent of L. This is illustrated by the two largest systems L = 128 
(diamonds) and 256 (solid line), with data just on top of each other. 

From the relations between T(k) and ^'(k) the above result is of relevance 
for g(k). We therefore conclude that our results for the correlation function 
are without significant finite size effects only if T ^ 0, which means that 
precisely this condition has to be fulfilled to facilitate reliable determinations 
of the correlation length ^. 

In the following determination of Tc and u we only make use of the data 
for the larger lattices, L = 128 and 256. For L = 128 the criterion T fti 
suggests making use of data for T/J > 0.84, only. The corresponding 
temperature limit for L = 256 is more difficult to obtain since T*^^^^) suffers 
from large statistical errors. The temperature limit we have used, T > 
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Figure 15: The correlation function g{q) for two different lattice sizes, together 
with the helicity modulus for the smaller size. Panel (a) is the q = component 
for system sizes L = 64 and 128 whereas panel (b) is the q = 27r/128 component 
for L = 128 and 256. In both cases g{q) start to differ when T 7^ 0. 
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Figure 16: Size dependence of T(k) at T/J = 0.840 for k = {kx,Q). The message 
in this figure is that T(k) and thereby g(k) is size-dependent unless T(k = 0) = 0. 



0.828 J, is obtained from considering the fraction L/\. The temperature 
hmit for L = 128 gives (L/A)min > 27, which for L = 256 is fulfilled only for 
T/J> 0.828. The determination of A is discussed in Sec. ^A . 

One should, of course, keep in mind that the usual finite size effect, related 
to the fraction L/C,, also may be relevant in these systems. It does, however, 
seem that this condition is the more restrictive one only for sizes L > 256. 



4.3 The correlation length in the FFXY model 

We now apply the methods and results from the previous Sections to our 
MC data for the FFXY model. 

Figure |l3(a) shows our values for ^ obtained by self-consistently fitting 
^(q) with q < tt/C and q < a/oT to Eq. (p]) for systems of size L = 64, 
128, and 256. The data for the larger sizes are also given in Table ^. Also 
shown is the helicity modulus for the two smaller sizes. We note that ^ 
for L = 64 (solid dots) start to deviate from the results for the larger lattice 
(L = 128, open squares) at the temperature where T^^^^ becomes appreciably 
different from zero. The corresponding situation holds for obtained with 
L = 128 and 256. Note also that this finite size effect is somewhat peculiar 
since analysis of data from a smaller lattice yields a too large value of the 
correlation length. This is opposite to the usual case, cf. Fig. |T^. 

After skipping the data affected by finite size effects, the analysis of the 
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correlation function in the FFXY model, becomes very similar to the corre- 
sponding analysis of the 2D Ising model in Sec. The data to be used is 
shown in Fig. p!7|(b). 

For determinations of Tc we first assume u = 1 and fit our data to Eq. 
([23| ) for temperatures T < Tmax- Our results for the critical temperature 
vs. Tinax; ^re showu as open circles in Fig. |l^(a). The figure shows a slowly 
increasing trend in for decreasing Tmax? similar to the results for the Ising 
model, Fig. |l^(a). From this we get our best value of the Z2 temperature, 
Te/J ^ 0.8225(5). The line in Fig. 0(b) is from the fit with T^^JJ = 0.842. 

The next step is to do a similar fit with z/ as a free parameter, by fitting 
to Eq. (|2^). Our values of and u are shown by solid squares as functions 



of Tmax in Fig. |TB|. For large T^ax the analysis gives non-Ising exponents, 
u ~ 0.9, but as Tinax decreases (the temperature interval shrinks) the data 
suggest an increasing trend in u towards 1.0. From this analysis of the 
temperature dependence of ^ we therefore conclude that the data, indeed, is 
consistent with u = 1. 

We also note that this dependence of u on Tmax is very similar to the 
corresponding analysis in the 2D Ising case. Fig. |l5(b). Analyses based on 
data from a somewhat too large temperature interval suggested u < 1.0, 
but with decreasing temperature interval the correct value of the exponent 
was obtained. This comparison serves to strenghten our conclusion of the 
ordinary 2D Ising value for the correlation length exponent in the FFXY 
model. 



4.4 Screening length A 



It is also of great interest to determine the screening length A, associated 
with the free vortices in the system. This turns out to, in some respects, be 
similar to the above determination of ^. The starting point is the expected 



behavior for the wave- vector dependent helicity modulus p2 



T(k) = 



J 



J 



e(k) ek2 + A-2^ 



where A is the screening length associated with free vortices and e is due to 
the polarization of bound pairs |]22|. Neglecting the fc-dependence in e one 



expects k^/T(A;) to be linear in k^. This is, however, not quite the case. 
The data shows a negative curvature and we therefore perform a fit to the 
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• L = 64 

*L = 128 
+ L = 256 




0.84 0.85 0.86 0.87 
T/J 

Figure 17: Correlation length ^ in the FFXy model, (a) Examination of the 
finite size effects. The data corroborates the conclusion that the determinations 
of ^ suffer from finite size effects if T 7^ 0. (b) The correlation length data to be 
used in the determinations of Tc and v. The solid line is obtained by fitting to Eq. 
(23) for T < 0.842J. Estimates of the statistical uncertainities in ^ are found in 
Table S. 
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Figure 18: Determinations of Tc and u from the temperature dependence of ^. 



This is a careful analysis of the data shown in Fig. 17(b) and listed in Table ^. 
The fits are performed with several different temperature intervals, including data 
up to Tmax- The open circles are obtained by fitting to Eq. (pSj), i.e. assuming that 
V = \. The solid squares are obtained with as a free parameter by fitting to Eq. 
(24). Panel (a) shows the obtained critical temperatures whereas panel (b) show 
the corresponding values of v. This is very similar to the behavior in the 2D Ising 



model, shown in Fig. 14. 
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Table 3: Correlation length versus temperature for system sizes L = 128 and 
256, together with rough estimates of the associated statistical errors. NMCS 
is the number of Monte Carlo sweeps through the system for the respective 
sizes. The data above the horizontal lines suffer from finite size effects and 



are therefore not used in the determinations of and v in Sec. 4.3 
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second-order polynomial 



k2 ^ ^ 

— flo + «ik^ + a2(k^Y. 



Tik) 



Assuming that the above equation holds for wavelengths larger than A we, 
self-consistently, make use of data for k < ir/X, only. The screening length is 



then obtained from A = yai/ao- The results are shown in Fig. |T^(a). Note, 
again that the finite size effects set in for L = 64 and 128, at temperatures 
below ~ 0.87 and ~ 0.84, respectively. 

We also note that the data obtained does fit well to the well-known 



Kosterlitz expression for the characteristic length 1 23 1. This is shown in Fig 



|T9|(b) where we plot In A vs. l/yT/TxT — 1- From this linear curve it is 
possible to get a value for the screening length at T^. An extrapolation to 
Tc/J = 0.8225 gives A(Tc) ~ 17.7. From our rough criterion for negligable 
finite size effects in Sec. |4.3| , L/X > 27, we may then estimate that the finite 
size effects associated with the KT transition would be negligable for systems 
with L > 478. 

The slope C ~ 0.42 is not quite in accordance with C ~ 0.54 obtained 
from the finite-size analysis of T/, in Sec. |3.4|. This kind of difference was 



also found in the ordinary XY model with no frustration! 27 1. However, it 
is only the analysis in Sec. |3.4| that probes the region immediately around 
Tkt, and this appears therefore to be the more reliable one when it comes 
to determining the asymptotic behavior. 



4.5 Different boundary conditions 

Figure |20| show g{r) versus r at T/J = 0.81, and 0.82, respectively, obtained 
with both PRC's and FBC's for lattice sizes L = 32, 64, and 128. At the lower 
temperature the results from larger lattices are squeezed in between the PBC 
and FBC results for a smaller system. When this behavior is valid, it seems 
safe to conclude that the behavior of an infinite system is somewhere between 
these two limits, and, to an excellent approximation, may be obtained as the 
average of these two curves. 

The behavior is dramatically different at the higher temperature. This 
is the case even though this temperature is well below Tc. (We expect this 
different behavior to set in for T > Tkt-) At this higher temperature, g{r) 
for both PRC's and FBC's decrease with increasing lattice size, and it is 
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Figure 19: Temperature dependence of the screening length, A. (a) The inverse 
screening length versus temperature. The finite size effect gives to large values of 
A. Considering the values from the largest lattices (the uppermost values) the data 
has a clear curvature. This is in contrast to the corresponding behavior of 1/^ in 



Fig. 17. This is a clear demonstration that A behaves very differently from ^. (b) 
Verification of the Kosterlitz temperature dependence for A. Here Tkt/<^ = 0.8108 



from Sec. 3.4 and the data affected by finite size effects is removed. The slope from 
this figure is, however, not in accordance with C = 0.54 it 0.02 from the finite-size 



scaling analysis in Sec. 3.4. 
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therefore not possible to obtain any safe estimate for the behavior in the 
thermodynamic hmit. 

We are now in the position to determine Moo(Tkt)- This is of great 
interest since our argument for two distinct transitions in section |3.3.2| was 
based on a simple way to estimate this quantity. Our values for Moo(T) are 
obtained by taking the average between data from PBC's and FBC's. M^q 



obtained in this way is shown in Fig. |2l|. With Tj^^/J = 0.8108 we find 
Moo(Tkt) ~ 0.783(2), which, indeed, is a good candidate to the L ^ oo 
limit in Fig. |. The results from our more elaborated analyses are thus in 
very good agreement with the simple approach of Sec. p.3.2| . 

4.6 The correlation length ^ for T < 

In this section we focus on the Z2 correlation length in the low-temperature 
region, T < Tc. As discussed above we actually need g{q) for an infinite 
system for a reliable determination of C,- As shown in Fig. pn|(a) the finite size 
effects may be virtually eliminated by taking the average of g{r) for PBC's 
and FBC's. The same turns out to be true for the Fourier components g{q)- 
However, as discussed above this only works at T < Tkt, which means that 
it does not seem possible to get any values for ^ right below Tc. 

The starting point for the present analysis is MC data obtained with both 
PBC's and FBC's for a lattice of size L = 128, in a temperature interval 
0.770 < T/J < 0.811. Much as in Sec. [4.3| we fit l/g{q) to an expansion in 



q^, cf. Eq. (pll). The data should be taken at small g; we restrict the analysis 
to qq < VO.l. The main difference compared to the high-temperature case, 
is that the data for g = has to be excluded at low temperatures. This is 
so since g{q = 0) is directly related to the staggered magnetization squared 
and g{q) therefore is not a smooth function at g = 0, in the low temperature 
phase. 

Figure ^ shows the obtained values of ^. There are two things to note. 
First, the data fall to a good approximation on a straight line, 1/^ oc (T— T*), 
but with the temperature T* significantly different from Tc ~ 0.8225J. The 
reason for this is possibly that Tkt lies between this data and Tc, and the 
system undergoes a dramatic change at the KT transition. 

The second point of interest is the value of the correlation length at Tkt- 
The finite size scaling analysis of T in Sec. p.4| would not be reliable if there 
were significant finite size effects associated with the Z2 degrees of freedom. 
We therefore need data with L ^ ^. From Fig. ^ we find ^(Tkt) ~ 7.2 
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Figure 20: The correlation function g{r) with different boundary conditions. 
Panel (a) is for T = 0.81J « Tkt, whereas panel (b) is for Tkt <T = 0.82J < T^.. 
At the lower temperature the correlation functions obtained with the two different 
boundary conditions scale with the system size in opposite ways, which facilitates 
a determination of the behavior in the thermodynamic limit. At the higher tem- 
perature this is no longer true and there is no easy way to extrapolate to the 
behavior in the thermodynamic limit. 
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Figure 21: The staggered magnetization. The figure shows the average of values 
obtained with PBC's and FBC's for L = 128. Since these sets of values are very 
close (difference < 0.01) we expect this average to be an excellent approximation 
of the behavior in the thermodynamic limit. 




T/J 

Figure 22: The Z2 correlation length ^ in the temperature region below Tc (and 
Trt)- Because of the dramatic change in the system at the KT transition, the 
critical temperature may not be determined from this kind of data. 
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which means that the systems used in the finite size scahng analysis, indeed, 
are considerably larger than the correlation length; L/C, > 4:A for the smallest 
system, L = 32. 



5 Discussion 



Our conclusion for the FFXY models with Villain or cosines interaction is 
thus that there are two distinct transitions, with Tkt < Tc. This is in contrast 
to the conclusions from a study of frustrated systems with a variable bond 
strength —7] J at one bond per plaquette||l8|. Here rj = 1 corresponds to the 



fully frustrated case, and both smaller and larger values of r] were used in 
order to separate the transitions, with the hope that this information would 
shed light on the behavior of the FFXy model, t] = 1. The conclusions - Fig. 
7 in Ref. [0 - were that as r] increases above rjc = 1/3 there are two distinct 
transitions - with < Tkt ^ which merge (or nearly so) at ?7 = 1. As 
is increased further the transitions separate again, still with the Ising-like 
transition at the lower temperature. 

Since the result at r] = 1 were not precise enough to exclude two distinct 
transitions, it is perfectly possible to harmonize these results with our present 
finding (Tkt < Tc for 77 = 1). The resulting picture is then the following. As 
a function of rj the temperatures Tkt and Tc approach each other, cross at rj 
below but close to 1, cross again above but close to r] = 1 and then separate 
as 7] increases. While this might be an unexpected scenario, it seems to be 
consistent with the conclusions based on the XY Ising model in Ref. |p4| . 

To summarize the findings of the present paper, we have found ample 
evidence for two distinct transitions with Tkt < Tc in the Villain version of 
the FFXy model. Our results corroborate the conclusions presented in Ref. 



T7| . We have given a strong argument for two distinct transitions based on 
the universal jump condition alone, which amounts to demonstrating that 
the staggered magnetization is finite, and actually quite large, at the KT 
transition. Since this result is obtained with no complicated analysis or 
detailed fitting whatsoever, we consider this to be a very robust argument 
for the existence of two transitions in the FFXF model. 

A detailed finite size scaling analysis of the helicity modulus gives at hand 
that the model undergoes an ordinary KT transition at Tkt/ J = 0.8108(1), 
and also that the size- dependence of this quantity in narrow regions both 
below and above Tkt is just as expected for a KT transition. Furthermore, 
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with this value for Tkt it is found that the temperature dependence of the 
screening length is in accordance with the well-known Kosterlitz' expression. 

From thorough studies of the correlation functions in both the FFXy 
model and the 2D Ising model, we have determined ^ and the correlation 
length exponent u. Our studies give at hand that the obtained values of u 
depend on the size of the temperature intervals from which ^ are taken, and 
that u approaches 1.0 as the size of the temperature interval is decreased. 
The result is therefore clearly in favor of ordinary Ising exponents in the 
FFXY model. Using u = 1 the Z2 critical temperature was determined to 
7;/^ = 0.8225(5). 
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A Duality transformation for a frustrated model 

The exact duality transformation [pO|] of the XY model with Villain interac- 
tion is between the XY model with PBC's and the Coulomb gas with a term 
proportional to the polarization squared pl|, |ll|. The derivation and nota- 
tions below closely follows Ref. fill. For an arbitrary set {Aij} the partition 
function is 



Z 



fU^ ( E m - - AM . (25) 



With the Fourier expansion of the Boltzmann factor 



-h^/2l3J 



the partition function becomes 

z= exp i ^ h,jA,j 1/ n ^ exp ^ n 



{hij}=-oo \ (ij) I I \ je{jl) I (ij) 



v/2^ 



Integration over 9i gives restrictions on the set {hij}, J2j hji = for each /. 
A new set of variables, defined at the centers of the plaquettes, take care of 
these restrictions: 

K = 'S'r+y/2 — 'S'r-y/2, 

h^r = Sr-±/2 - 5'r+x/2, (26) 

together with S = {S^, Sy) which are included to allow for J2r K = 
y ||1 1[| . To take care of the frustration we make use of D x Ar = 27r/r. 

^ h^A^, + h^Al = E('5r+y/2 - 5r-y/2)A^ + (^r-x/2 - 5r+x/2)A^ 
r r 

= X! 'S'r(^r-y/2 ~ ^r+y/2 + ^r+x/2 ^ ^r-x/2) 
r 

= 5:5r(Dx Ar)=27r5]5r/r. 
r r 

/ \ -[5. -5,,] 2/2/3 J 

52 



using 



this becomes 



CI rm 



cr 



(r,r') 



2I3J 



{mr},m 



i2Ti m ■ (J + ^(/r + mr)crr 



where m = (mx, ^^j,)- After a Fourier transform this may be written 



W (c/RecTk) ((ilmcrk)(icr exp 

k>0 



y2/3JL2 ^ G'(k) ~ 



1 



^ exp 

{mr},m 



L2 



(27) 



where M = m + X^r "^rf is the polarisation. After a transformation back 
to ordinary space, this gives 



Z = exp 

{mr},m 



1 TVT^ 
I - 2 EK + /r)G(r - r')K. + U>) + — 



Dropping the 'spin wave part' Zgw, speciahzing to the frustration /r = 1/2, 
and substituting mr + /r ^ ""^d we may then write 



{mr=±l/2...},m 



1 



(28) 

which is our final result for the partition function. 

We then turn to demonstrating the equivalence of Eqs. (^ and (|^). The 
wave- vector dependent helicity modulus is defined through p5[ 



T(k) 



dAldA'_^ 



where the derivative is with respect to an that is transverse, i.e. k ■ A^ = 
0. Since only the rotation of A contributes to the change in free energy, we 
make use of 2-7? /k = ik x Ak to obtain 



T(k) 



k2 „ d^F 



47r2 9/k9/_k 
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From Eq. (P7| ) then follows 
1 d^Z lA-K^J 



|m„k + /-k|G(k)|mk + /k|G'(-k) 



which, with F = —TlnZ, = — 1/G'(k), and rrik + /k ^ m^, reproduces 
Eq. (i. 

By instead starting from Eq. (|25|) we find 



1 d^Z 
Zdf^dU 



dAr 



df-v 



To evaluate the derivatives we have to choose Ar such that D x Ar = 27r/r. 
One choice that fulfills this is 



271 ~ ~ 

Ak = fk.{iky, —ikx 



that gives 



2n 



{iky, ikxje 



dh k2L2 

Putting all this together and introducing the current jr = {U'{(f)f),U'{(l)^), 
the wave-vector dependent helicity modulus becomes 



and with 27rJfr = D x jr -^^ 2-71 Jv^ = ik x jk, we finally obtain Eq. (^. 
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